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1. Introduction 

Cosmic rays (CRs) are charged particles of extraterrestrial origin impinging on the Earth’s 
atmosphere. They are mostly made up of protons, with a small fraction of other fully ionized 
nuclei, electrons, and antiparticles. The CR energy spectrum is smooth and almost featureless from 
a few GeVto several EeV. It can be described by a power law with two breaks: the knee (steepening 
at a few PeV) and the ankle (flattening at a few EeV). 

CRs are deflected by (extra-)galactic magnetic fields, making it impossible to identify their 
sources (acceleration sites) directly for the energies considered here. However, we can identify 
acceleration sites by studying neutral byproducts of CR acceleration, such as y-rays and neutrinos. 
Another way to study the sources of cosmic rays is to precisely measure their spectrum and com¬ 
position, and compare these to the predictions from models of different types of acceleration sites, 
see e.g. [1]. Examples of such candidate acceleration sites are supernova remnants for galactic 
cosmic rays and gamma-ray bursts or active galactic nuclei for extra-galactic cosmic rays. Some 
of the results of these precision measurements are shown in Eigure 1 . 

1.1 Experimental Methods in Cosmic Ray Physics 

In the MeV to TeV range, cosmic rays can be measured directly by balloon- or satellite-borne 
detectors placed at the top of the atmosphere, see for example [2]. The energy, mass and charge of 
the cosmic ray can be determined directly, from the energy deposited in the detector and the shape 
of the particle’s track. 

Upon interaction with the Earth’s atmosphere, cosmic rays with energies above a few hundred 
TeV will produce a cascade or shower of secondary particles (nucleons, electrons, positrons, pions, 
muons, photons, neutrinos). If the primary particle’s energy is high enough, part of this extensive 
air shower (EAS) may be detected by arrays of particle detectors on the ground. Energy and direc¬ 
tion of the primary particle can be inferred from the density and arrival times of the secondaries. 

Eor intermediate energy primaries, no, or only a few, secondary particles will reach the ground. 
However, the charged component of the shower will emit Cherenkov radiation in the visible/UV 
range, high up in the atmosphere. This light can be detected by Imaging Atmospheric Cherenkov 
Telescopes (lACTs), mostly used for y-ray astronomy. These instruments use very sensitive cam¬ 
eras to image air showers and reconstruct energy, direction, and type of the primary particle from 
the brightness, shape, and orientation of the image. Using the direct Cherenkov technique [4] 
described below, they can measure the charge of the primary particles, especially for heavy nuclei. 
This enables the use of large existing data sets for cosmic ray physics (typically up to 1000 hours of 
quality-selected data per year of operation, with a few hundred cosmic ray events per second, which 
are considered as background for y-ray astronomy). The systematic uncertainties are dominated by 
the atmosphere and thus largely complementary to the direct detection and EAS experiments. 

1.2 Direct Cherenkov Technique 

In addition to the Cherenkov light radiated by the charged component of the shower, charged 
primaries with high velocities will also radiate direct Cherenkov (DC) light before starting a shower. 
The DC light is emitted at high altitudes under very small angles. Eor today’s generation of lACTs, 
with typical pixel sizes on the order of 0.15° [5], this light is generally concentrated within one (or 
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Figure 1: Energy spectra of various elements in cosmic rays, measured by different experiments over a large 
energy range, from [3]. 


two neighboring) pixel(s). The intensity is proportional to Z^, the squared charge of the primary 
nucleus. In combination with an energy estimate from the intensity of the shower image, the DC 
light can be used to obtain separate spectra for light and heavy nuclei if the DC contribution can 
be identified as such and distinguished from the light emitted by the secondary particles. Light 
elements such as helium or oxygen are not very well suited for the direct Cherenkov technique 
as the DC light emitted by these elements is too dim compared to the contribution from the air 
shower. For example, iron emits ten times as much DC light as oxygen (Z = 8). In the following, 
iron (Z = 26) will be used as a representative for heavy elements as it is the most abundant element 
in cosmic rays with Z > 20. 

The DC technique is sensitive to iron nuclei in cosmic rays with energies of about 10 TeV 
to several 100 TeV. For lower energies, the particle will interact and start a shower before emit¬ 
ting Cherenkov light. At higher energies, the light emitted by the shower is brighter than the DC 
contribution. For more details, see [4; 6; 7]. 

1.3 The VERITAS Experiment 

VERITAS (Very Energetic Radiation Imaging Telescope Array System) [5; 8] is an array of 
four lACTs located at the Fred Lawrence Whipple Observatory (FLWO) in southern Arizona (31 
40N, 110 57W, 1.3km a.s.L). It has been operational since 2007. Each telescope has a total mirror 
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area of ~100 m^, and a camera consisting of 499 pixels (photo multiplier tubes), each with a field 
of view of 0.15° diameter. It is sensitive to showers induced by y-rays from 85 GeV to > 30 TeV. 

The standard data analysis chain is described in [9, 10]. After charge integration and image 
cleaning, a moments analysis is performed and the Hillasparameters [11] are calculated. From the 
orientation and shape of the images, the direction of the primary particles and the position of the 
shower core on the ground can be determined, and y-ray induced showers separated from cosmic- 
ray induced ones. Look-up tables are used to determine the energy given the intensity of the shower 
and the distance to the telescope. These results may be used directly to determine spectra, fluxes, 
etc, or used as input to more advanced reconstruction techniques such as the template likelihood 
fitting method described below. 

2. Template Likelihood analysis method 

This reconstruction method requires a model that predicts the average photon intensity in the 
camera, depending on some properties of the primary particle such as energy, arrival direction, 
etc. The probability distribution of the signal (charge) per camera pixel is then determined by 
those parameters as well as the detector response and its uncertainty. Given the images taken of 
a particular shower, this probability distribution can be converted into a likelihood function of the 
primary parameters. Maximizing this likelihood function produces a set of parameters which are 
good estimates for the “true” properties of the primary particle. In addition, the goodness-of-fit for 
those parameters can be used to separate background events (that are not described by the model) 
from signal events (that are well-described by the model). 

A similar likelihood fitting method has been used for the reconstruction of y-ray induced show¬ 
ers by other experiments, see for example [12], as well as by VERITAS [13], to improve the 
performance of the experiment. 

In this paper, we describe the performance of a variant of the likelihood fitting method used 
to reconstruct the energy of cosmic ray irons from their associated air showers, and to separate 
iron-induced air showers from those induced by lighter cosmic rays. The model parameters are 
the primary particle’s energy E, arrival direction in camera-centered coordinates (A^, T^), height of 
first interaction h and position of the shower core on the ground {Xp, Yp). The light intensity in the 
camera is modeled as a set of templates, using Monte Carlo simulations of iron-induced showers in 
a fixed grid in E, h, and the distance between the telescope and the shower core. Several simulated 
showers are averaged to obtain a stable average of the photon intensity in the camera. Interpolation 
is used to predict the light intensity for arbitrary parameter values. Rotation and translation of the 
template image account for different arrival directions and core positions on the ground. 

CORSIKA [14] was used for air-shower simulation including the emission of Cherenkov light, 
and the grisudet package^ was used for ray-tracing in the telescopes. Figure 2 shows an example of 
the predicted average light intensity in the camera for an iron shower. The contributions from both 
DC light and the air shower can be clearly seen. A more detailed description of the reconstruction 
process, can be found in [13]. The likelihood function adapted for the reconstruction of iron 
showers is described in [15]. 

'http://www.physics.Utah.edu/gammaray/ GrISU/ 
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VERITAS —ICRC 2015 
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Figure 2: Predicted average light distribution in the camera {template) for iron showers with E — 30 TeV, 
with a distance of 50 m between the detector and shower core, h = 33 km. The origin is defined by the 
primary particle’s arrival direction. The x-axis is dehned by the long image axis. Note the contribution from 
direct Cherenkov light at around 0.1° from the primary’s direction, clearly distinguishable from the broader 
contribution from the Cherenkov light emitted by the shower particles. 


3. Method Performance 

To test the method, a set of templates was produced using simulated iron showers from zenith 
in the energy range from 10 - 100 TeV. These templates were tested against a set of simulated 
iron showers with random energies, core positions and arrival directions in the camera, which was 
passed through a simulation of the full VERITAS detector, again using the grisudet package. These 
showers were analyzed with the template fitting method described earlier. Starting values for the 
template fit were obtained using the standard geometrical reconstruction, with special lookup tables 
for iron-induced showers. 

The following three quantities are plotted in Figures 3-8 and compared: The energy bias, 
defined as fhe relafive error in fhe energy — 1, fhe energy resolufion (68% confainmenf infer- 
val around fhe median reconsfrucfed energy), and fhe angular resolution (68% confainmenf radius 
of fhe reconsfrucfed direcfion around fhe frue direcfion). Energy bias and resolufion are imporfanf 
sources of sysfemafic uncerfainfies on fhe energy specfrum. The angular resolution is ploffed here fo 
show fhaf fhe direction is reconstructed well enough to identify the DC light contribution. All quan¬ 
tities are plotted both for the template likelihood fitting reconstruction (“likelihood”, solid lines) 
and the standard reconstruction (“standard”, dashed lines). The template method always performs 
better than the geometric reconstruction. Quality cuts are applied, but no cut on the goodness- 
of-fit after the likelihood minimization. Note that the performance of the standard reconstruction 
is worse for iron-induced showers than for y-ray induced showers (as described in [8]). This is 
due to the fact that hadronic showers tend to be much less smooth, with larger shower-to-shower 
variations than y-ray induced showers. 

Figure 3 shows the energy bias depending on the arrival direction (offset from the camera 
center), with the telescopes pointing at zenith. Above 50 TeV, it the energy bias is very flat. There 
is a significant dependence on the offset from the camera center: For large offsets, the bias is as 
large as 10%. This happens because the templates were produced for zero offset. The VERITAS 
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optical point-spread function degrades with the offset, so there is less light and hence the energy is 
underestimated. Since the energy bias is relatively flat, this is can be corrected for in the future. 

Figures 4 and 5 show the energy resolution and angular resolution for different offsets. The 
resolution is significantly better than for the geometric reconstruction, and the dependence of the 
energy resolution on the offset it greatly reduced. However, below about 80 TeV, the angular 
resolution is still larger than the pixel diameter, which implies that the DC light, which should be 
contained within one pixel, is not always identified properly. 

Figures 6, 7, and 8 show the corresponding plots for the telescopes pointing 20° away from 
zenith, averaged over all azimuth directions. The energy bias is slightly worse, due to the fact 
that the light had to travel longer through the atmosphere, so more light was absorbed. Energy 
and angular resolution are not degraded. The dependence on the azimuth was investigated as well 
(not shown). There is no significant dependence of the energy bias, energy resolution or angular 
resolution on the azimuth. 

4. Conclusions and Outlook 

To fully understand the sources and acceleration mechanisms of cosmic rays, it is important to 
measure their composition and spectrum over a large energy range. The range of tens to hundreds 
of TeV is not well-covered by existing direct detection experiments. However, lACTs are sensi¬ 
tive to cosmic rays in that energy range. By detecting direct Cherenkov light emitted high in the 
atmosphere, they can separate light and heavy cosmic rays and produce separate spectra. 

A method for reconstructing the energy of iron-induced showers, imaged by the VERITAS 
experiment, has been presented. This likelihood fitting technique relies on fitting template images 
to recorded data, and can be used to separate iron-induced showers from background events as 
well. However, further work is needed to improve the performance of the fit and to optimize the 
separation between showers induced by iron and those induced by light nuclei. 

In the future, this method will be applied to archival data sets taken with the VERITAS in¬ 
strument. The goal is to obtain an energy spectrum of iron nuclei in cosmic rays for energies in 
the range of about 30 to 300 TeV. Due to the improved reconstruction method, this is expected to 
reduce the uncertainty on the spectral measurements compared to current results. 
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Figure 3: Energy bias for 0° zenith angle. 
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Figure 6 : Energy bias for 20° zenith angle. 
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Figure 4: Energy resolution for 0° zenith angle. 


Figure 7: Energy resolution for 20° zenith angle. 
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Figure 5: Angular resolution for 0° zenith angle. 


Figure 8: Angular resolution for 20° zenith angle. 
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